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In this paper we apply Kraichnan's direct interaction approximation, which is a one loop pertur- 
bation expansion, to magnetohydrodynamic turbulence. By substituting the energy spectra both 
from kolmogorov-like MHD turbulence phenomenology and a generalization of Dobrowolny et al.'s 
model we obtain Kolmogorov's and Kraichnan's constant for MHD turbulence. We find that the 
constants depend of the Alfven ratio and normalized cross helicity; the dependence has been studied 
here. We also demonstrate the inverse cascade of magnetic energy for Kolmogorov-like models. Our 
results are in general agreement with the earlier simulation results except for large normalized cross 
helicity. 
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I. INTRODUCTION 



Turbulence still remains primarily an unsolved problem. Among existing statistical theories of fluid turbulence, 
Kolmogorov's phenomenology is the most prominent one. Kolmogorov |5J hypothesized that the energy spectrum 
E(k) of fluid turbulence in the inertial range, i.e., for scales between the energy feeding and dissipation range, is given 
by 

E(k) = K Ko U 2 / 3 k- 5 / 3 , (1) 

where Kko is an universal constant called Kolmogorov's constant, k is the wavenumber, and II is the nonlinear cascade 
of energy, also equal to the dissipation rate of the fluid. This power law has been confirmed by experiments Q and 
numerical results || . 

Later on various theories were developed to understand fluid turbulence, the primary ones being Direct Interaction 
Approximation (DIA) of Kraichnan 0,(5 1, Wyld's field-theoretic technique [[|, renormalization groups (RG) [Q, mode- 
coupling Eddy-Damped Quasi Normal Markovian (EDQNM) closure schemes ||. All these theories yield results 
consistent with the Kolmogorov's power law. These calculations yield K ss 1.5 in three dimensions (3D). In two- 



dimensions (2D) Kraichnan 10 1, Olla |ll|], and Nandy and Bhattacharjee [[t2| show that Kko ~ 6.4 in region where 
inverse cascade of energy occurs. 

The DIA of Kraichnan is one of the fully consistent analytical turbulence theory of fluid turbulence. It 
essentially involves perturbation theory similar to that used in quantum field theory. In this calculation Green's 
function and correlation function are calculated self consistently to first nonvanishing order. Only triad interactions 
are allowed to this order; for this reason it is called direct interaction approximation (refer to Figure 1). The 
DIA of Kraichnan Q yielded k~ 3 ^ 2 energy spectrum for fluid turbulence, which was found to inconsistent with the 
experiments. However, later on it was shown that DIA can also yield k~ 5 ^ 3 energy spectrum ||. The RG and mode 
coupling theories are similar to DIA, and they have been shown to be yield results close to those from DIA. 

For magnetohydrodynamic (MHD) turbulence also there are several phenomenologies. In MHD we have velocity 
field u and magnetic field B = Bo + b, where Bo denotes the mean magnetic field and b denotes the magnetic field 
fluctuations. In this paper, in place of u and b, we use Elsasser variables z ± = u ± h/y/inp, where p is the density of 
the plasma. The wave speed of Alfven waves due to the mean magnetic field, called the Alfven speed, is denoted by 
Ca and is equal to Bo/-\/47rp. The variables z^ denote fluctuations having positive and negative velocity-magnetic 
field correlations respectively. Throughout the paper we assume that p is constant, which implies that the fluid is 
incompressible. Marsch, Matthaeus and Zhou, and Zhou and Matthaeus |l3| proposed a phenomenology similar to 
Kolmogorov's fluid turbulence phenomenology. We call this Komogorov-like MHD turbulence phenomenology. In this 
phenomenology the energy spectra E (k) of the z fluctuations are 

E ± (k) = k± (n+) 4/3 (n-)~ 2/3 k- b ' 3 (2) 

where are the nonlinear energy cascades of fluctuations, and K are constants, whom we refer to as 
Kolmogorov's constants for MHD turbulence. The spectral normalized cross helicity a c (k), defined as (E + (k) — 
E~ (k)) I \E + (k) +E~(k)), and spectral Alfven ratio r^(/c), defined as E u /E b , play an important role in MHD turbu- 
lence. Here E u and E b are the kinetic and magnetic energies respectively. In this paper a c (k) and r^(fc) are taken to 
be independent of k, and we will usually drop the term spectral while referring to these constants. According to the 
assumptions of the model, this phenomenology is expected to be applicable when amplitudes of the fluctuations are 
greater than the mean magnetic field or the magnetic field of largest eddies, i.e., z 3> i?o- In this paper we theoret- 
ically calculate the values of K . The quantity E R = [E u — i? b ), called the residual energy, plays an important role 
in MHD turbulence, as will be seen later in this paper. 

Kraichnan [jl4| argued that when the mean magnetic field or the magnetic field of the largest eddies Ca is much 
larger than the fluctuations z , the magnetic energy spectrum E b (k) is given by 

E\k) = A(U b C A ) 1/2 k- 3/2 , (3) 

where n h is the magnetic energy cascade rate, and A is Kraichnan's constant. Dobrowolny et al. jl5| generalized this 
model for z ± and found that 

n+ = n - = A^C^E+i^E-i^k 3 , (4) 

where A is a constant. In this paper we also calculate the value of the constant A for a modified Dobrowolny et al.'s 
model. 
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Matthaeus and Zhou, and Zhou and Matthaeus developed a generalization that contains both Kolmogorov-like 
and Dobrowolny et al.'s models as limiting cases. In their model the dissipation rates are given by 

n± = A-*B^E+(k)E-(k)k* 
C A + ^kEHk) 

In this paper we do not discuss this model because of its complexity. 

Some of the analytical theories of fluid turbulence have been applied to MHD turbulence. For example, Pouquct 
et al. and Grappin et al. [[[7) applied the EDQNM closure schemes, and Fournier et al. and Camargo and Tasso |l(f | 
have applied the RG technique to MHD turbulence. Verma and Bhattacharjee |l8| have applied DIA technique to 
MHD turbulence when E~(k) = E + (k) and obtained the Kolmogorov's constant K . In this paper we have generalized 
their formalism for cases when E~(k) =^ E + (k) and obtained K . In this paper we have also theoretically obtained 
Kraichnan's constant A for a modified Dobrowolny et al.'s model. 

The solar wind is an ideal natural experiment where some of the predictions of MHD turbulence theories could be 
tested. The solar wind observations by the spacecraft show that energy spectra of the solar wind is closer to /c~ 5 / 3 
rather than to fc~ 3 / 2 jn|. This is in spite of the fact that the mean magnetic field in the solar wind (approximately 5 
nanotesla) is approximately 2 to 10 times larger than the large-scale fluctuations. However, it should be noted that 
the exponents 5/3 and 3/2 are difficult to distinguish. Recently Verma et al. [^0) have performed direct numerical 
simulations of MHD turbulence in which they could not conclude whether the spectral index was 5/3 or 3/2. However, 
in their simulation the nonlinear energy cascade rates 11^ appear to follow Eq. (|J) rather than Eq. (H|). Note that 
the EDQNM calculations of Grappin et al. [[l7| show that the spectral index varies between and 3 depending on 
normalized cross helicity. 

The outline of the paper is as follows: In section 2 of this paper we derive the equations for the correlations functions 
and flux functions that are used in subsequent sections. In section 3 we show that the energy spectra of Eq. (||) are 
consistent solutions of MHD equations and find the constants i"T . We find the constants for various E~(k)/E + (k) 
and E M (k)/E b (k) ratios. In section 4 we construct a generalization of Dobrowolny et al.'s model, and then calculate 
the Kraichnan's constant that model. Section 5 contains discussion and conclusions. Other than a brief remarks in 
sections 4, we do not attempt to address the important question under what ratio of z^ /Bq the transitions from the 
energy spectra fc -5 / 3 to fc~ 3 / 2 take place; the answer to this question require further analysis which is beyond the 
scope of this paper. 



II. MHD TURBULENCE 



The MHD equations in terms of are |1 



^-z ± T C A ■ Vz ± + z=F • Vz ± = --Vp tot + ^+V 2 z ± + z/_V 2 z T + f ± , (6) 
at p 



iy±c 2 /(Ana) 

v± = 2 (?) 

where C A = Bo/V47rp, p tot is the total pressure, i.e., thermal plus magnetic pressure, f ± are forcing, v is the kinematic 
viscosity, c is the speed of light, and a is conductivity. The corresponding equation in Fourier space is 

(i + v+k 2 ) zf(k, t) + v-Pz?(k, t) T i{C A ■ k)z±(k, t) = 

/±(k) - ik tP tot (k, t)-J dp Pj zf(p, t)z±(k - p, t) W 

We can eliminate pressure p tot using the incompressibility condition, which is 

hzf(k,t) = 0. 

After the elimination of pressure the incompressible equation is 

(i + v+k 2 ) z±(k, t) + v-Pzf{k, t) T i{C A ■ k)z±(k, t) = f . 

/±(k) - eM ijm (k) J dpzf(p, i)4(k - p, t) W 

where 
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M ijm (k) = kjPirnik); P im (k) = S im - (10) 

We have introduced an expansion parameter e which is set to one at the end. We assume that the forcing / + and 
/~ are independent, i.e., df Sl /df S2 = S SlS2 , where si, S2 — ±. Note that when z + = z = u, we recover the Navier- 
Stokes equation. The quantity Mjj m (k) for MHD is different from that in fluid turbulence because of asymmetry 
introduced by z + and z . 

Now we define z — z correlation tensors Sij as follows: 

<z+(k,t)z+(k',0) = S++(k,t,t')S(k + k'), (ii) 
(zr( k ,i)z-(k',0) = Sy-(k,M')<5(k + k'), (12) 

(zf(k,t)zr(k',t')) = S+-(k,t,t')S(k + k'). (13) 

Here (x) denotes an ensemble average of x. We assume that the turbulence is stationary, which implies that the 
correlations arc independent of time: S^j S2 (k,t,t) = S^ S2 (k). For fluid turbulence Orszag |)| derived correlation 
functions and their properties. Our following discussion is similar to that of Orszag. 
From Eq. (0) we can obtain equations for correlation tensors that are 



(JL+2V+&) S± ± (k,t,t) + 2v_k 2 S+-(k,t,t) - 2i(C A ■ k)S±±(k,M) 
2F±(k) - eM Um (k) J dpT| m (-k, p) - eM, 7m (-k) / dpT± n (k, p) 



(14) 



where 



(£ t +2v + k 2 )S+-(k,t,t) + v_k 2 (S++(k,t,t) +ST j -(k,t,t)) 

-2i(C A -k)S+-(k,t,t) = 2F*(k) (15) 
-eM iim (k) J dpT™(-k, p) - eM jlm (-k) J dpT« (k, p) 



(z±(k,i)zf(p,i)z±(q,i)) = T± m (k,p)<5(k + p + q ) (16) 

(« i -(k,t)z J -(p,t>+(q,t))=2^(k,p)*(k + p + q) (17) 

(4(k,^+( Pl t>-(q,t)> = Z$»(k,p)*(k + p + q) (18) 

</±(k)z±(k'))=F ; ± m (k)J(k + k') (19) 
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- (/+(k)z J "(k')+/r(k>+(k')) = F^(k)«5(k + k') (20) 



For isotropic turbulence 



S s i r(k)=P l3 (k)C s ^(k). (21) 



In presence of a mean magnetic field, the correlation tensors are expected to be anisotropic. In the solar wind 
the fluctuations are anisotropic; the amplitudes of the fluctuations along the Parker field is one-third that of those 
perpendicular to the Parker field pi]] . However, the analysis of anisotropic turbulence is quite complicated, and in 
this paper we assume isotropy for all the correlations. 
We also define one-dimensional energy spectra as 
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2v + k 2 )E ± (k,t,t)+2v_k 2 E R (k,t,t)=T ± (k,t) (24) 



For isotropic distribution in k space, after substitution of i = j in Eqs. ( [i^ , [T5| ) and taking the trace, we obtain 

'd_ 
dt 

2v + k 2 ) E R (k,t,t) +z/_fe 2 (E+(k,t,t) +E-(k,t,t)) = T R {k,t) (25) 



d 
dt 



where 



T±(fc,t) = -4irk 2 M t3m (k)Im J dpT± m (k,p) (26) 
T R (k,t) = -^k 2 M tjm (k)Im /dp ( ^( k .P)-^(k, P ) \ _ (27) 



Note that the term (C^ • k)E ± (k) = C Ak s E ± {k) 1 where k s is the wavevector along the mean magnetic field, is zero 
for isotropic turbulence because it is an odd function of k s . 

When v+ = V- = 0, one can easily check from the above equations that the total energy = J E ± (k)dk is 
conserved. However, total residual energy E = J E R (k)dk is not conserved. In terms of T's these statements of 
energy conservation are 



oo 

± 



T ± (k,t)dk = (28) 



o 



(29) 



but 

/ T R (k,t)dk^Q. 
Jo 

We can also derive a theorem of "detailed conservation of energy E ± (k)" , i.e., 

s ± (k,p,q) + s ± (k,q,p) +s ± (p,q,k) + s ± (p,k,q) + s ± (q,k,p) +s ± (q,p,k) = 0, (30) 

where 

s ± (k, p, q)5(k + p + q) = -7m ((k ■ Z T(p)) (z±(k) ■ z ± (q)) ) (31) 

The diagrammatic representation of the terms S + (-, •, •) is shown in Figure 1. The solid lines represent z + and the 
wavy line represent z~ . For S~(-, •, •) we need to change the solid lines to wavy lines and vice versa. The physical 
interpretation of these diagrams is that the energy from any z + (k) mode gets transferred to modes z + (p) and z~(q) 
(k + q = k), and so does energy from z~(k). However, the quantities E + (k)+ E + (p)+ E + (q) and E~(k)+ E~(p)+ 
E~(q) are separately conserved. 

In presence of viscosity there is dissipation in the plasma. The dissipation rates of z fluctuations are given by 



e 



± _ 



-f t J™E±{k,t)dk 



= 2v + f °° k 2 E ± (k, t)dk + 2v_ J °° k 2 E R (k, t)dk. 



(32) 



Integrating both sides of Eqs. (|24|,|25|) over a sphere of radius K in k space, we obtain the energy cascade rates 
H±> R (K,t) that is 

n±(x, t) = - J K r±(fc, t)dk = /~ r±(fc, t)dk 

f t J K E ± {k,t)dk-2^ + J^°k 2 E ± {k,t)dk-2^J^k 2 E ± {k,t)dk. 1 ' 



Il R (K,t)=-f*T R (k,t)dk 
= -Aj^ E R (k,t)dk-2is + f™k 2 E ± (k,t)dk-2v + J™k 2 E ± (k,t)dk. K ' 

The above equations will be used in the subsequent sections. In Eqs. (|2jj2^) we find that the nth-order moments 
involves (n + l)th-order moments along with nth and lower order moments. Hence one cannot close these equations. 
This is the famous "closure problem" 91. To circumvent this difficulty, various techniques have been devised, primary 
ones being renormalization groups DIA mode-coupling scheme ||, EDQNM closure schemes |^Jl^] etc. 

In this paper we only discuss DIA of MHD turbulence. 
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III. DIA OF MHD WITH Ca = 



In this section we assume that the mean magnetic field or the magnetic field due to the large eddies is small as 
compared to the fluctuations, i.e., Ca <C zM. Here we set Ca = 0. In the next section_we will analyse for the case 
when Ca 3> zM. Firstly, we calculate the Green's function in the spirit of Kraichnan 



A. Computation of Green's Function 

It is convenient to derive Green's function in Fourier space k = (k, u>). Following Leslie Q we define Green's function 



(35) 



In Fourier space the equation of Green's function is 



-iuj + u + k 2 v^k 2 
v_k 2 —iuj + v + k 2 



dp 



p+q=k 



G+ + (k) G±-(k) 

z i (p) z+(q) 



= P ii (k)-cM iim (k)5 



Gi 3 + {q) Gj-rtf) 



(36) 



where p = (p,o/) and q = (k — p, u>— u>'). We denote 



— iuj + v + k 2 
v_k 2 



—iu> + v + k 2 



= G°(k, w ) 



(37) 



Now we solve the Eqs. 



G in (k,w) = P in (k)G(fc) 
and (pq) perturbatively. We expand and Green's function G in series: 



.MO) ft 



(k) 



ez 



±(i) 



(k) + 0(e 2 



(38) 



(39) 



G£*»(fc) = G°r W (k) + <r [1 '(k) + 0(e 2 ) 



SlS 2 (l), 



(40) 



In DIA we keep terms only up to first nonvanishing order J5|. By substituting these terms in Eqs. (|^) and (|36|), and 
equating the terms with equal powers of e, we obtain 



*r (0) (*) 



= G°(k) 



J % 

fr 



(41) 



\k) 



-iG u {k)M tjk (k)x / dp 

Jp+q=k 



zf 0) (p)z m ^(q) 



-(0), 



(42) 



G (1) (k) 



-G u (k)M ljm (k) 



dp 



p+q=k 



{a \p) z+ (0) 

+ (0) 



(<?) 



++(o) 

mi 

-+(o) 



(q) g+7 (0) (9) 
(p) G;r {0) (p) 



(43) 



The factor 1/2 in the equation appears because of trace of Py(fc). In DIA f is assumed to have a gaussian distribution. 
Therefore, according to Eq. (f4l|), zM°) also have a Gaussian distribution. Green's function G to leading order is 
(refer to Eq. |||) 



G 



G(k) = G<°> (k) - §e 2 4 0) (k)M ljm (k) / _ dp 

Jp+q=k 



Z 7 (0) (p) ^ (0) (<7) 
Z 7 (0) (p) ^ (0) «) 



(44) 



In the above equation we have used the fact that G^(k, w) is statistically sharp ||, and 
(Z^G^) + 0{e 2 ). The substitution of G 1 of Eq. @ in the above equation yields 



,±(0) 



0, (ZG) 



G(k) = G°(k)--G° il (k)M ljm (k) 



dp < 



p+q=k 



L z7 (0) (p) z-r(q) 



-(o), 



A/ sa6 (<?)G+t (0) (-z) / A/ sab (g)G+; (0) (g) / <//• 

J f+s— g ./ f+s— q 

M safc (p)G7+ (0) (p) / df Af safc (p)G7 s - (0) (p) [ 



J r-\-s—p 

z ; (0 \f)p bn (s) z+ (Q \s)p an (f) 



(0), 



r+s— p 

G++(°)(s) G+7(°)(s) 



<V (0) (0 G7 r (°)(f) 



> 



(45) 



Following Kraichnan's DIA we substitute e — 1 and replace z^ by z and G^ * 1 in 2nd and 3rd line by G. 
Comparing the resulting equation with 



G(k) = G°(k) - G°(k)t(k)G°(k) 



where £ is self-energy, we obtain 

E(fe) = 



£-+(&) s — (fe) 



V +(k) V-(k) 

n+(k) r,-(k) 



(46) 



(47) 



r] + (k) = \k 2 J fi+g= ' k dpx 
[h(k,p,q) (G-+(p) + G—(p)) C+-(q) + b 2 (k,p,q) (G++(p) + G+-(f>)) G" (<?)] 



(48) 



[b 3 (k,p, q) (G-+(p) + G—(p)) C++(q) + h(k,p, q) (C++ (p) + G+"(p)) C+~ (q)] 
Here bi(k,p,q) — k~ 2 Bi(k,p,q), and B^s are ]j| 



Bi{k 7 p,q) 
B 2 (k 1 p,q) 
B 3 {k,p,q) 
B 4 {k,p,q) 



kp(—z + z 3 + x 2 z + xyz 2 ) 
kq(l + z 2 )(z + xy) 

kp(—Z + z 3 + xy + x 2 z + y 2 z + xyz 2 ) 
kp(—z + z 2 + y 2 z + xyz 2 ) 



(49) 



where (x, y, z) are the cosines of angles between (p, q), (q, k), and (k, p) respectively. 

When z+ = z~, we have G ++ = G~~ = G H , G ++ = G~~, and G H = G h . Under this condition we expect 
S ++ = £~~, but Eqs. ( p7| , p8| , p9| ) do not yield this results. However, the equality can be easily achieved by sym- 
metrizing Eqs. (h_8|49j) by interchange of + and — signs in Eq. (45) and averaging of the resulting equations with 
corresponding Eqs. (|48|) and (|4S|). The final equations for rf+ are 



r]+(k) = lim^o \k 2 J dpx 
X [b x {k,p, q) (G++( P ) + G++(p)) C+-(q) 
+b 2 (k,p,q) (G++(p) + G++(p)) C++(q) 
+b 3 (k,p, q) (C++ (p) + G++ (p) ) C++ (q) 
+b A {k, Pl q) (G*=F(p) + G^±(p)) G+"(g)] 



(50) 



Now when z+ = z , we recover the fluid limit. In this case C++ = G~ = C+~ and £++ = S~~ = 2T, fluld . 
Here wc take a; — > limit in our calculation; this approximation is valid for the large and inertial scale fluctuations. 
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After this we solve for E and G self consistently. A typical self-energy digram is shown in Figure 2. Considering that 
z + = z correspond to fluid limit, we postulate the following relaxation behaviour for the z — z correlation functions: 

C^% U ) = ^-^W ;C *F(fc,a,) = -L °y , (51) 



C±±(fc, - t x ) = C±(fc) exp (-4 (t a - tO) ; 
C±T(fc,t 2 - tl ) = C*(fc) exp (-^i (t 2 - ti)) ; 



The Eq. (|5J) yields 



(52) 



G(fc) 



(53) 



which after substitution in Eq. (|50| ) and integration of a/ integral @ yields 



r ] ± {k)=jdp 



•A; 2 



Vp+Vp + Vq+Vq 

b 3 (k,p,q)C^ (q)r,± 
+Vp +2i)J 



fe 4 (fc,P,g)C^(g)^ 

Vp 1 +Vp +Vt +Vq 



(54) 



We have ignored the contributions from the pole at u>' — because it corresponds to a static situation. The quantities 
77 + /2 and rf~ /2 correspond to the inverse of the response time of the eddies z + and z~ respectively. 

Our self-energy matrix is differs from that obtained by Camargo and Tasso's using renormalization group calcula- 
tions Jl(|; their £ is symmetrical with both the diagonal elements equal. This is due to different assumptions of the 
models. In Camargo and Tasso's model cross helicity is zero which forces both the diagonal elements to be equal. 
However, our model can be applied to situations where cross helicity is not zero. Also, we assume independent / + and 
/~ forcing, whereas Camargo and Tasso assume independent f u and f h forcing. Renormalization group calculation 
for varying normalized cross helicity is under progress. 

We need to evaluate the integrals in the above equation In 3D the integral || is 



dp 



dpdq 



2npq 



(55) 



while in 2D it is 



dp 



dpdq 



sin(p, q) 



The integral is to be performed in the hatched region of Figure 3. 

We substitute the power-law expression of C + (k),C~(k), and C R (k) from Eqs. 
interpreted as the inverse of time-scale of z ± , we substitute B 



V 



: (k) = A ± (rF) 2/3 (n 



r ±\-V3 



(56) 



g3j23j). For 17* (fc), which is 



(57) 



where A ± are constants. We also assume that the inertial range spectra of the residual energy E R (k) is proportional 
to A: -5 / 3 . The EDQNM calculations of Pouquet et al. and Grappin et al. jl7| show that E R {k) is proportional to 
k~ 2 . Marsch and Tu analysed the solar wind data and found that the residual energy spectra is between /c~ 3//2 
and £;~ 5 / 3 , but closer to k~ 3 ^ 2 . Hence, there is no consistent theory for the spectra of the residual energy. In our 
model, for consistency, we assume that E R (k) oc fc~ 5 / 3 . With this power law the ratio between E R (k) and the total 
energy E (k) is constant; we denote this constant by a. Similarly the ratio E~ (k) / E + (fc) is a constant in our model; 
we denote it by [3. 



Now the substitution of power laws for rfj: in Eq. ( pi} ) yields Q 



(a±t 



fc 2 (l^,K) + &3(l,^,K) 
,(l+C ±1 )? 2 / 3 + 2C± 1 K 2 /3 



J J d^dnqn 



-8/3, 



2(1+£T1) 



&i(l,<,,K)+b 4 (i,<;,«) 
(i+C ±1 )(<; 2 / 3 +k 2 / 3 ) 



(58) 
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in three dimensions. Here £ = (A~II + ) / (A+Il - ) ,p — C,k and q = nk. In 2D the terms uw~ 8 / 3 and wv~ s ^ 3 are 
replaced by u> -8 / 3 and v~ 8 ^ 3 respectively, and the whole integrand is divided by 7r(l — x 2 ) 1 ^ 2 [^|Jl0[|. 

The integrals of Eq. (E^) suffer from the well known "infrared problem" which comes from the strong dynamic 
coupling of fluctuations with widely differing wavenumbers. Over the years, various techniques have been developed 
to tackle the difficulty in the context of pure fluid turbulence. Simplest among all these methods is introduction 
of cutoff, which is discussed in Leslie [||. Later on Renormalization group technique 0], the Lagrangian or semi- 
Lagrangian pictures p2| ] , and self-consistent screening ]8| were used to make the integral naturally finite and cut off 
independent. The infrarred difficulties associated with Eq. ( |58| ) can be similarly resolved. Knowing that the full 
theory is constrained to be finite, we adopt the practical procedure of evaluating the integral in Eq. ( f47| ) with a cut off 
ko = Xk and choosing A — 1 so that the pure fluid value of A 2 /K are obtained correctly when a = (3 = 1. Thereafter 
A is not varied. 

In 2D when E u < E b , we choose A = 1, same as A of 3D, but not 0.065 which yields K — 6.6 in fluid limit. The 
choice of same A in 2D MHD was motivated by fact that in MHD turbulence, forward cascade of energy and inverse 
cascade of cross helicity occur in both 2D and 3D (refer to Ting et al. p3|]). Note that in 2D fluid turbulence an 
inverse cascade of energy occurs, hence the behaviour of fluids turbulence in 2D and in 3D are dramatically different. 
It appears that A = 1 is not applicable for cases when E u > E b ; therefore we have calculated the constants for 
E u > E b . 

To obtain the numerical value of if and A^ 1 we need two additional equation involving and A^ . The additional 
equations we use are the equations of cascade rates of z^; these equations are derived in the next subsection. 



B. Computation of Cascade rates 

The cascade rates or the fluxes have already been defined in terms of T ±,R respectively (ref. Eqs. d33"| , |34|) ) . 

We evaluate T ±,R to first order, and then, following Kraichnan |||| substitute G°, C°'s by G's and G's respectively 
(refer to Leslie [0 for this procedure applied in fluid turbulence). We will symbolically illustrate the evaluation of 
T up to first order. For this part of the calculation we work in (k, t) space (refer to Leslie ( 0) for details). The 
quantity T is proportional to equal time triple correlations. To first order T is 

T(k, t, t) <x / p+q=k dp (z(°>(k,t)z( )(p,i)z«(q, t)) 
+ (z( )(k,i)z( 1 )(p,i)z( )(q,i)) + (z( 1 )(k,i)z( )(p,i)z(°)(q,<)). ^ 

Note that (z^(k,t)z^(p,t)z^(c[,t)) is zero. Now using Eq. ( [f2"| ) we expand z^ and substitute in the above 
equation, which yields 

T(k,t,t) oc J k dp l^drf^dt'x 
<z( )(k,t)^ (0) (P^)G( 0) (q^-^)^ (0) (i-,i')2 (0) (s^')) + •■• 



T(k,t,t)oc dp dt'G (0) (q,t-t')G(k,t-t')G(p,t-t') + --- (61) 

«/p+q— k J — oo 

The terms G's and C's appear as G SlS2 and C SlS2 in the final expression. The fnll expression for is 

T±(k,t,t) = 4nk*f p+ci=k dpf t _ oo dt>x 
[-6i(A,p, q) {G^(p, At)C R (q, At)C±(p, At) + G^(p, At)C R {k, At)G±(q, At)} 
-b 3 (k,p, q) {G^ip, At)C R (k, At)C ± (g, At) + G^(p, At)C R (q, Ai)C ± (fc, At) J 
-b 2 (k, q,p) {G^iq, At)C R ( P , At)C R (k, At) + G^lq, At)C^(p, At)G±(fc, At)} (62) 
-b 4 {k, q,p) \G ±z f (q, At)C*(p, AtfC+ik, At) + G ±± (g, At)C R (p, At)C R (k, At)} 
+b 2 (k, q, P ) (G ± T(fc, At)C R (p, At)C R {q, At) + G ±± (A:, At)C T (p, At)G ± (g, At)} 
+h(k,p, q) {G±T(fc, At)G T (p, At)G±(g, At) + G±±(fc, At)C R ( P , At)C R (q, At)}} 

where At — t — t f . The expression for T R is similar to T . Substitution of T's yield the following cascade rates: 

p OC p P 

U^ R (k)= dk' dpdqS^ik'.p^q) (63) 
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where the dpdq integral is over the hatched region of Figure 3. Note however that the hatched region A has q > p. 
If we have p > k, then all three vectors k, p,and q forming triangle will have magnitudes greater than k; hence these 
triad will not contribute to the flux. Therefore, in region A the flux contribution comes from the following region of 
integration: 



/>oo />k rp+k' 

U ±M (k)= dk' dp dgS ±lK (fc' ) p,g) 

J k JO J p* 



where 



(64) 



For region B similar arguments lead to 

ir^jfe) 



p* = max(p', \p — k'\). 



i rk rq+k' 

dk' dq 

JO J q* 



dpS^ H {k',p,q) 



where q* = max(g', \ q — k'\). 

After substitution of C ±,R and G in the Eq. (61) and with the change of variable || 



we obtain 



. k vk wk 
k = -;p = — ;q = — 

u u u 



A± 



K+K- 



dv In ■ 



l+v 



dw 



dw In ■ 



l+tu 



dv 



* ± (l,u,w) 



(65) 



(66) 



(67) 



(68) 



where in 3D 



ttA _ l n+K+R- 

11 - 2 AT 



f dv In - f „ + " dw + f dw In — f 1 ~J~ 1 " dv 



I rl+w 



* ± (1,W, w) 



(1+ g T l) [bl{l,V,W) +b 3 (l,V, W)) X 



_ t 1 ^ , 1 r 

(2+(l + 5±!)(u 2 /3 +t0 2/3^ ^ ((l+C±!)(l+i; 2 /3) + 2^2/3) 

~f r ^ 1 (b 2 (l,w,v)+b,(l,w,v))x 

1 , a 2 (l+/3) 2 /(4/3) 



(2+2C± 1 t; 2 / 3 + (l+C± 1 )u) 2 / 3 ) T ((l+C± 1 )(l+« 2 / 3 +t« 2 / 3 )) 

+^ff^(b 1 (l,v,w) + b 2 (l,w,v))x 

_ 1 ^ , a 2 (l+/3) 2 /(4/3) 

(l+e± 1 +2C± 1 i> 2 / 3 +2«) 2 / 3 ) ~ ((l+5± 1 )(l+-u 2 /3+ t « 2 / 3 )) 



(69) 



(70) 



-8/3 



a 2 (l+/3) 2 /(4/3) 



rrry (6i(l,U, w) + 6 3 (1, V,W)) X 



((l+0(l+t> 2 / 3 +tu 2 / 3 )) (2C+(l+C)i) 2 / 3 +2tu 2 / 3 ) 

- "" 8/ Yi+o +/3)/2 (62(1, + 64(1, «;,«)) x 
+ 



(2f+2£-u 2 / 3 + (l+£)tu 2 / 3 ) ((l+?)(l+«> 2/3 )+2^ 2 / 3 ) 

+ ( ( "i+ c -i / ) 3 v,w) + 62(1, I0,u)) X 

1 , a 2 (l+/3) 2 /(4/3) 

(l+C+25t' 2 / 3 +2«; 2 / 3 ) ((l+C)(l+i; 2 / 3 +io 2 / 3 )) 



(71) 
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* W2 (1, v,w) = (bi(l,v,w) + h(l,v,w)) 

g 2 (l+/3) 2 /(4/3) 



((l+C)(l+i> 2 / 3 +™ 2 / 3 )) (2+(l+4)v 2/3 +2£u> 2 / 3 ) 

- " (1+^-1)2 f ( & 2(l,w,i;) + 6 4 (1, tu,«)) 



(2+(l+£)(« 2/3 +™ 2 / 3 )) ((l+?)(l+M! 2 / 3 )+2{t> 2 / 3 ) 

+ ( 7i+e) 7 ( fo l(M,w) + fo 2(l,W,D)) 

' (l+^) 2 /(4;3) 



(72) 



(l+£+2i> 2 / 3 +2£w; 2 / 3 ) ((l+C)(l+w 2 / 3 +io 2 / 3 )) 

In 2D the terms ww" 8 '' 3 and wv~ 8 ^ 3 are replaced by w -8 / 3 and v -8 / 3 respectively, and the whole integrand is divided 
by tt(1 - a; 2 ) 1 / 2 [M. 



C. Computation of K and A* 

We numerically solve for if ± and A 1 * 1 from Eqs. (Bq, |6^) for given ratios of a = E R (k)/E(k) and /3 = E~(k)/E + (k). 
Since the integrand itself is a function of f, which it self is a function of A 1 * 1 ,]!" 1 ", and II - , we start with a guessed 
value of £ and iterate until the two successive values of £ are approximately equal. We perform these calculations for 
both 2D and 3D. The values of K ± , A ± , and £ for various a and (3 for 3D are listed in Table 1, whereas for 2D, they 
are listed in Table 2. 



D. Computation of K R 

If we postulate that 

E R (k) = K m (n R ) 2/ V 5 / 3 , 

then the substitution of H R in Eq. ([39]) yields 

E R (k) \ f A~ 



2/3 



(73) 



(74) 



where I is the value of the integral of Eq. ( |69| ) . The values of K Rl computed for various a and (3 are listed in Tables 
1 and 2. 

We can model E R (k) in yet another way. This modeling is motivated by Eq. (|l5|). From this equation we argue 
that 



Therefore, 



or, 



£ fl (/c) 



£*(fe)fe 



2 ^JO fc -6/3 n Jl 



v ^k+ (it-) 2/3 (n+r 1/3 + Vi^(n+) 2/3 (n-r 1/3 ' 

The constants i4T fl2 calculated using the above equation are listed in Tables 1 and 2. 



(75) 



(76) 



(77) 
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E. Discussion 



Several important conclusions can be drawn from the numbers shown in Tables 1 and 2. The Kolmogorov's constant 
for MHD turbulence are not universal constants as compared to universal fluid Kolmogorov's constant. We find 
that the constants depend on Alfven ratio ta and ratio of E~(k) and E + (k). They could also depend on the 
mean magnetic field, but study of this effect is beyond the scope of this paper. 

As reported in our earlier paper Jl8|], when (3 — E~(k)/E + (k) = 1, the constants K + = K~ = K increases 
monotonically from 1.43 to ~4.07 as we go from fully fluid case ta = oo to ta ~ 0.2, then it decreases and finally 
reaches 3.51 when the plasma is fully magnetic. Please note that the numbers here are slightly different from those 
in JL8[| ; the discrepancy is due to the lower accuracy used in the earlier paper. In 2D when (3 = 1, K is in the range 
of 5 — 7.5 for < < 1. In both 2D and 3D for a given ta, as we decrease (3 (or increase normalized cross helicity 
cr c ), K + (constant corresponding to larger of E + and E~) decreases whereas K~ increases. 

The above DIA results are in qualitative agreement with the preliminary numerical results of Verma et al. [po[ . 
Verma et al. performed direct numerical simulation of MHD turbulence in 2D (8 runs) and 3D (single run) for various 
initial normalized cross helicity and mean magnetic field. When E + ps E~ , the constant K + was approximately equal 
to K~ . The mean values of the constant obtained by them were K = 3.7 in three-dimensions for E~ /E + = 0.6, and 
K = 6.6 in two-dimensions for E~ /E + = 1, 0.6. However, the constant K of the majority species (larger of z ) was 
always marginally lower than that of minority species. Unfortunately, we do not have the inertial range Alfven ratio 
TA{k) for our simulations; we will study the ta dependence in future simulations. However, we supsect the inertial 
range Alfven ratio to be in the range of 0.25 — 1. We find that for small a c the mean values of K from the simulations 
are generally close to the DIA result for ta — 1 and a c = 0. 

However, for large a c there appears to be serious discrepancies between the simulation results and our results. In 
simulations for large a c , n+/n~ ~ E~(k)/E + (k), however, in our DIA calculations n + /Tl~ <C E + (k) / E~ (k) in fluid 
dominated cases and n + /n~ ^ E + (k) / E~ (k) in equipartioned (ta — 1) or magnetically dominated case. Regarding 
constants K ± , for a c (k) ps 0.9(/3 = 0.05) Verma et al. found that the average K + ps 2.4 and K~ ps 24.0. In DIA 
calculations, for magnetically dominated cases, the equations for do not yield real solutions. This indicates that 
our model is not fully consistent atleast in magnetically dominated limit; this point is illustrated in the following 
paragraph. Nevertheless, there are some quantitative agreement between the simulation results and DIA results, e.g., 
the constant K corresponding to the majority species is always greater than the K corresponding to the minority 
species. 

When we turn off the velocity field, i. e., z + = z = b, the MHD equations are linear. Here the magnetic 
energy dissipates only due to the resistivity, and there is no turbulent cascade rate. However, in our formalism, we 
obtain nonzero turbulent cascade rate for this case. Also, there is no consistent solution for when a c is large and 
ta < I- To circumvent this inconsistencies we may have to modify our relaxation time, Green's functions etc., and 
also possibly may need to use u and b rather z . This work is under progress. 

Note that the ratio I\ R /\] + increases monotonically as r a decreases from infinity to zero. By the time r a = 5, the 
ratio n fl /n + has already become large. This large ratio implies that n b pa — II" (H u > and n b < 0), i.e., the kinetic 
energy is cascading from small k to large k, while the magnetic energy is cascading from large k to small k. This is 
a theoretical demonstration of inverse cascade of magnetic energy. Earlier work on amplification of magnetic energy 
have been done by Fyfe et al. using numerical simulations, Leorat et al. using EDQNM closure scheme, Gloguen 
et al. using scalar model of MHD turbulence, and Ting et al. by turbulent relaxation ( [^3| and references therein). 
Our results may be useful for modeling dynamos from turbulent energy cascade point of view. If we start with small 
magnetic field in a fluid dominated plasma, there will be a small flux of magnetic energy cascading to large-scales 
from small-scales; that will enhance the magnetic field in the system and decrease ta- As ta decreases, the flux 
cascading inversely also increases, hence further enhancing the magnetic field of the system. After the magnetic field 
has become sufficiently large {ta = 5), almost all the fluid energy flux goes into magnetic energy flux, which makes 
the enhancement of the magnetic field faster. A quantitative model based on these ideas could be useful in estimating 
the time-scales etc. in dynamos. 

In the above discussion we have not differentiated between inertial range ta (E u (k) / E b (k)) and the ta = E u /E b , 
which may be a gross approximation. Also modeling of relaxation of C H (fc, At) as well as modeling of E R (k) is not 
on a strong footing. We do not get consistent answers for equipartioned and magnetically dominated cases. Numerical 
simulations might provide important clues that could help us modeling these quantities better. 

In the next section we will apply DIA to MHD turbulence in situations when Ca ^S> z ± . 
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IV. DIA OF MHD WITH C A > Z 



As mentioned in the introduction, in presence of a strong mean magnetic field, the energy spectra is expected to 
be proportional to k~ 3 l 2 . However, the solar wind observations |l9| as well as the numerical simulations [^0) do not 
appear to support this hypothesis at least for Ca/z ± < 5. It is possible that when Ca/z ± 3> 1, e.g., when the ratio 
is of the order of 100, the Kraichnan's |L4| or Dobrowolny et al.'s (l5) arguments are valid. Unfortunately numerical 
simulations have not yet probed in these regimes. 

There exist another possibility that Bo OT C 'a which should be substituted in Eq. (^) is not the mean magnetic 
field but the renormalized Ca ■ When the mean magnetic field term of Eq. (||) is comparable to the nonlinear term, 
simple dimensional arguments yields the following scaling for Ca'- 

C A =Il 1/3 k~ 1/3 . (78) 

where II is the total energy cascade rate. This argument is similar to that due to Bhattacharjee |2j| for renormalization 
of sound speed in randomly stirred fluids. We conjecture that the small-scale fluctuations are affected by renomalized 
magnetic field rather than the mean magnetic field. Hence, we should substitute renormalized Ca of Eq. ( |78|) in Eq. 
(||) that immediately yields k~ 5 ^ 3 energy spectra for energy. Therefore, large Ca could yield k~ 5 ^ 3 energy spectra 
rather than fc -3 / 2 as argued by Kraichnan |l4| and Dobrowolny et al. jl5). However, in absence of a quantitative 
theory or definite numerical results, we will make certain assumptions and calculate the Kraichnan's constant for a 
variation of Dobrowolny et al.'s model. 

As argued by Kraichnan Jl4| and Dobrowolny et al. JTj| we assume that the Alfven time-scale, rather that the 
nonlinear time-scale, is the relaxation time-scale when Ca/z^ ^S> 1. Using random phase approximation Veltri et al. 
f24| showed that Green's function and relaxation of correlation functions are 



exp [-gkC A (t2 




h)} 

exp [-gkCA(t2 - h) 



(79) 



C^ H (k, t 2 - h) = C^ik, h = t 2 ) exp [-gkC A {t 2 - t x )] . (80) 

We assume g = 1. Unfortunately, Veltri et al. assumed equipartition of magnetic and fluid energy which is not 
the case in our model discussed in this section. In Dobrowolny et al.'s model |ll| the dissipation rates n + and n~ 
are equal irrespective of E~ (k) / E + (k) ratios. This has been found to be inconsistent with the numerical results [po[ . 
To bring in some effects of E~(k)/E + (k) we modify the Dobrowolny et al.'s model in the following manner: 

E ± (k) = A ± (U ± C A ) 1/2 k- 3 ^ 2 (81) 

We do not need to solve for self-energy in this case; Green's function has been obtained using random phase 
approximation. We, however, solve for cascade rates. The equation corresponding to Eq. ( |62| ) is 

T±(M,t) = 4 7 rfc 4 / p+q=k dp/: oo ^'x 
[-bi(k,p, At)C R (k 7 At)C±(q, At) 

-h(k,p, q)G(p, At)C R (q, At)C±(fe, At) 

-b 2 (k, q,p)G(q, At)C=F(p, At)C±(fc, At) (82) 
-b 4 (k, q,p)G(q, At)C R (p, At)C R (k, At) 
+b 2 (k, q,p)G(k, At)C*(p, *t)C±{q, At) 
+h(k,p, q)G{k, At)C R ( Pl At)C R {q, At)] 

where G(k,At) = exp(— CAt). From the Eqs. (|82| ) and a similar equation for T R we can derive expressions n ± and 
IL R . 

From the equations for and H R we can solve for A using manipulations similar to those used in the previous 
section. The values of for various values of a and are listed in Tables 1 and 2. When (3=1, the constants 
A + = A~ = A are approximately 2.0-2.5 in 3D and 3.5-4.0 in 2D. The constant A obtained here for 3D is close to 



the one obtained by Matthaeus and Zhou |13|. As we vary j3, the constants A + and A~ begin to differ. We find that 
when E~(k) < E + (k), A + is always larger than A~ contrary to the Kolmogorov-like models in which K + < K~ . For 
ta < 1, the ratio n + /n~ is less than one, whereas for ta > 1, the ratio is greater than one. Hence ta appears to 
have significant effect on cascade rates in this model as well. 

The ratio n^/n -1 " decreases from a very large value (2519) to nearly zero as ta varies from zero to one. This implies 
that there is an inverse cascade of magnetic energy for ta < 1- As ta increases from 1, the ratio n fi /n + decreases 
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from zero to negative values (the magnitude of the ratio increases) till ta reaches around 15, then it again starts 
increasing. It reaches nearly zero at ta — 1000 and reaches one at ta = oo. This is an indication of inverse cascade of 
kinetic energy for 1 < ta < 1000. Qualitatively, when the magnetic energy dominates the plasma, there is an inverse 
cascade of magnetic energy. On the contrary, when the kinetic energy dominates (only till ta — 1000), there is an 
inverse cascade of kinetic energy. These predictions differ from those in section 3. Recall that in Kolmogorov-like 
models of the previous section, for all ta (except at oo) we had an inverse cascade of magnetic energy. 

Our generalized Kraichnan model is heuristic. Numerical simulations could help us in determining which of the 
MHD turbulence models are applicable in various parameter space. Comparisons of the numerical results with the 
above predictions will provide insights into these puzzles. 

V. CONCLUSIONS 

In this paper we have applied direct-interaction approximation to MHD turbulence and obtained the values of 
the Kolmogorov's constants for MHD and Kraichnan's constant for various E~(k)/E + (k) and E u (k)/E b (k) ratios. 
Ours is a field-theoretical (nonvanishing) first order calculation. We have obtained equations for the constants as 
well as ratios of energy cascade rates for both Kolmogorov's and Kraichnan's power laws. We have solved for these 
constants for various values of E~{k)/E + {k) and E u (k)/E b (k) ratios. When a c = 0, in three-dimensions the constants 
K + = K~ = K are close to 2 for fluid dominated cases, but between 3 and 4 for cases when ta < 1- In 2D, however, 
K is in the range of 5 — 7.5 for < ta < 1. For a given ta as we increase a c (or decrease j3), K + (the constant 
corresponding to majority species) decreases whereas K~ increases. 

The DIA results for Kolmogorov-like phenomenology with small a c are in good agreement with the simulation 
results of Verma et al. |2C|j . However, the DIA results for large a c are only qualitatively consistent with the numerical 
results. For example, the Kolmogorov's constant for (K + here) the majority species is always smaller than that for 
minority species, consistent with the numerical results. However, for large a c our DIA calculations yield n + /n~ <C 
E + (k)/E~(k) in fluid dominated cases and n + /Tl - ^> E + (k)/E~(k) in magnetically dominated cases; these results 
are in disagreement with the simulation results where n + /n~ ~ E + (k)/E~(k). We believe that the inconsistency 
can probably be removed if we modify the relaxation time or E R expressions, or possibly we may have to use u 
and b varibles. However, note that in 3D simulation only a single run with a relatively lower resolution (128 3 ) was 
performed. We need to perform more runs to compare the DIA results with the simulations. 

In this paper for the first time we have calculated Kraichnan's constants for a generalization of Dobrowolny et al.'s 
model. We find that when (3—1, the constant A + = A~ = A is approximately 2.0-2.5 in 3D and 3.5-4.0 in 2D, 
consistent with earlier estimates of A. As (3 decreases, the constants A + and A~ begin to differ. The DIA results for 
Kraichnan's model are inconsistent with the results of the numerical simulations performed so far. It is possible that 
Kraichnan's model or its generalizations becomes applicable for large Ca/z ± 1 however, at present we do not know the 
critical Ca/z where transitions from Kolmogorov-like models to Kraichnan's model take place. Another possibility 
is that C a of Kraichnan's model is the renormalized mean magnetic field which is scale dependent (oc fc -1 / 3 ), and 
that will yield fc -5 / 3 energy spectra. 

Our DIA calculations using Kolmogorov-like energy spectra exhibits an inverse cascade of magnetic energy flux. We 
also find that as r a decreases, the inverse cascade rate of magnetic energy increases. This is a theoretical demonstration 
of inverse cascade of magnetic energy, and this result could be useful for modeling dynamos from turbulent energy 
cascade point of view. However, DIA calculations using the generalized Dobrowolny et al.'s model yields an inverse 
cascade of magnetic energy for ta < 1, and an inverse cascade of kinetic energy for 1 < ta < 1000. 

The results discussed in this paper is useful for various applications in the solar wind. For example, with the 
knowledge of the Kolmogorov's constants for MHD, Verma et al. |2(| have calculated the contribution of turbulent 
dissipation to the overall heating in the solar wind. We anticipate that similar theoretical calculations can be useful 
in understanding of inverse cascade of energy, process called dynamo mechanism, and inverse cascade of cross helicity, 
a process called dynamic alignment. 
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Figure 1: Diagramatic representation of the nonlinear interaction terms s + 



(k-z+(a))(z-(b)-Vz+(c)), 



where a, b, c = (k, p, q) (refer to Eq. (30)). The solid lines represent z + and the wavy lines represent 



Figure 2: Diagramatic representation of self-energy S(fc,w). The quantities G and C in the figure are dressed 
quantities. 



Figure 3: The region of integration / _ k 
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TABLE I. For 3D MHD turbulence Kolmogorov's constants K , A , A etc. for various values of ta and E / E + . 
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TABLE II. For 2D MHD turbulence Kolmogorov's constants K ± ,A ± ,A ± etc. for various values of r A and E / E + . The 
constant A which appears in the lower limit of the self-energy integral is chosen to be 1.0. 
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